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SUMMARY 


In NASA Contractor Report 3^63 [Ref. l], maximum like- 
lihood and constrained least-squares methods were developed 
for estimating the parameters of the von Karman transverse 
and longitudinal spectra and autocorrelation functions. 
However, during the application of these methods to re- 
corded turbulence measurements [NASA Contractor Report 
3^64, Ref. 2] it was discovered that a small percentage of 
such records are better modeled by the Bullen spectra or 
autocorrelation functions. The present report extends the 
above mentioned maximum likelihood and constrained least- 
squares methods to Bullen transverse and longitudinal 
spectra and autocorrelation functions. 
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INTRODUCTION 


In Ref. 1, a general methodology was developed for 
characterizing atmospheric turbulence velocity records, 
estimating the parameters of the characterisations from 
such records, and computing aircraft response statistics 
from the characterizations. The characterization methods 
were implemented in Ref. 2 using several turbulence 
velocity records measured during the NASA MAT (Measurement 
of Atmospheric Turbulence) program. The turbulence model 
used in Refs. 1 and 2 assumes that a turbulence velocity 
record w(t) can be expressed as the sum of "slow" and 
’’fast" components w g (t) and w f (t) - i.e. 

w(t) * w s (t) + w f (t) 

= w g (t) + o f (t)z(t), (1.1) 


where 


w f (t) = o f (t)z(t) 

with 


( 1 . 2 ) 


a f (t) * 0 (1.3) 

and 


E{z( t ) } = 0, E{z 2 (t)} = 1, (l.U) 


where E{ • • * } denotes the mathematical expectation or en- 
semble average value of the quantity within the braces. 

The three random processes {w g (t)} , (o f (t)}, and (z(t)} 

are assumed to be mutually statistically independent. Fur- 
thermore, (z(t)} is assumed to be an ergodic Gaussian pro- 
cess. The process {o^(t)}is not assumed to be Gaussian. 

Realizations of (o f (t)} are assumed to vary slowly in com- 
parison with realizations of (s(t)}, and realizations of 

(w (t)} are assumed to vary slowly in comparison with 
s 


2 


ORIGINAL PAOt 14 
OF POOR QUALITY 


realizations of (w f (t)}. Hence, w Q (t) has the appearance of 
a slowly varying additive component to the process w f (t), and 
w f (t) has the appearance of ordinary turbulence with a tem- 
porally varying envelope that is governed by the behavior of 
Oj.(t). A more complete description of the above model can 

be found in Ref. 1. 

It has been shown in Ref. 3 that whenever the scale of 
fluctuations of o f (t) is of the order of ten or more times 

the integral scale of z(t), the velocity spectrum of w^(t) 

will, for practical purposes, coincide with the velocity 
spectrum of z(t) except for a constant multiplicative factor. 
Thus, if z(t) is regarded as ordinary turbulence with a 
von Karman transverse or longitudinal spectrum, then to 
within an excellent approximation Wj.(t) will possess this 

same spectral form. Hence, in Refs. 1 and 2, it has teen 
assumed that the velocity spectral densities or autocorrela- 
tion functions of w(t) are the superposition of the spectral 

densities or autocorrelation functions of w_(t) and 

. . s 

a* z(t), where o* denotes a temporal average or expectation. 

This assumption can be justified from the assumed mutual 
statistical independence of {w (t)}, (o f .(t)}, and (z(t)} 

and the slowly varying behavior of (o f (t)} relative to {z(t)}. 

Since the spectrum of w (t), by assumption, contains 

s 

predominantly low frequency or long wavenumber components, 
the asymptotic slope of the spectrum of w(t) should be gov- 
erned by the process z(t) which has been assumed to be de- 
scribed by the von Karman spectral form. The von Karman 
transverse and longitudinal spectral forms both possess the 
asymptotic -5/3 slope when plotted on log-log coordinates 
as required by classical turbulence theory [4], However, 
when the methods developed in Ref. 1 were applied to 7 real 
turbulence records in Ref. 2, it was found that the asymp- 
totic slopes of the spectra of 2 out of the 7 records were 
distinctly different from -5/3* These are the spectra of 
the lateral turbulence records from Flight 32 Run 4 and 
Flight 30 Run 8 of the MAT program shown in Figs. 5 and 11 
respectively in Ref. 2. Thus, an extension of the methods 
developed in Ref. 1 to accommodate asymptotic slopes dif- 
ferent from -5/3 would seem to be useful. 


3 
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Two methods were developed In Ref. 1 for estimating 

"spectrum" properties of velocity records. The first method 

is applicable to cases where the slow component w_(t) in 

8 

Eq. (1.1) is negligible in comparison with w f (t). This 

first method utilizes the maximum likelihood method to 
estimate the two parameters of the von Karman transverse or 
longitudinal spectral forms. The second method is appli- 
cable to cases where w_(t) in Eq. (1.1) is not negligible. 

s 

This second method itilizes a constrained least-squares inte- 
gral fit of an autocorrelation function model to the auto- 
correlation function of the measured turbulence record. 

Both of these methods are extended in this report to models 
where the asymptotic slope of the spectrum of the component 
z(t) is not constrained to be the classical value of -5/3. 


4 
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MAXIMUM LIKELIHOOn ESTIMATION OF PARAMETERS 
IN BULLEN TURBULENCE VELOCITY SPECTRA 

In Sec. 3 cf Ref. 1, the maximum likelihood method was 
i used to develop equations for estimating the integral scale 

and intens ty of the von Karman transverse and longitudinal 
spectral forms. In this section, this method is extended to 
the Bullen spectral forms given on pp. 200-202 of Ref. 5. 

The Bullen spectral forms are determined by three parameters 
rather than two as in the von Karman case. The third parame- 
ter controls the asymptotic slope of the spectra when plotted 
on log-leg coordinates. When this asymptotic slope is set 
equal to -5/3, the Bullen spectral forms reduce to the 
von Karman forms. 

The Bullen transverse and longitudinal two-sided wave- 
number spectra are [p. 200 of Ref. 5] 


and 


^ BT ^ ^ 


= o 2 L 


1+8tt 2 fc 2 k 2 (n+1 ) 
(l+4TT 2 S, 2 k 2 )^ + 3 / ’ 2 




a 2 L- 


(l+4'rr 2 fi. 2 k 2 ) n+1/2 


( 2 . 1 ) 


( 2 . 2 ) 


where the above forms are one-half of the values given in 
Ref. 5 because * m (k) and <t> T (k) are two-sided spectra rather 

D 1 t5L> 

than one-sided as in Ref. 5, and where the parameter a in 
Ref. 5 has been replaced by l to avoid later notational con- 
fusion. L is the scale of the turbulence, o 2 is the mean- 
square value, and £ is related to n and L by 


l 


T(n) L 

/fT r(n+Tp 


(2.3) 


where T ( a ) . is the gamma function. When n = 1/3, Eqs. (2.1) 
and (2.2) reduce to the von Karman forms given by Eqs. (3*18) 
and (3.19) of Ref. 1. 
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In the maximum likelihood procedure to be developed 
below, we shall derive likelihood equations for the three 
parameters a 2 L, l , and n In Eqs. (2.1) and (2.2). When 
solutions are arrived at for these three parameters, one 
can. readily compute from n and l 

r(n+i) 

L « £_ a , (2.14) 

T(n) 

according to Eq. (2.3), and then a 2 = a 2 L/L, which yield 
the scale L and mean-square value o 2 . 

Let Sj, j * 1,2,‘**,N denote the periodogram (unsmoothed 
spectrum) samples described on pp. 7^-80 of Ref. 1. The 
mathematical expectation of a typical spectrum sample S. 
can be expressed as J 


E(Sj) i: Sj = o 2 LFjU,n) , j=l,2,*‘*N (2.5) 

where for the Bullen transverse spectrum one has from Eq. (2.1) 


Fj U,n) 


1+ 8tt 2 2 k j (n+1 ) 

( i + i|TT 2 ?„ 2 k 2 ) n+ 3/ 2 

.1 


( 2 . 6 ) 


and for the Bullen longitudinal spectrum one has from Eq. (2.2) 


FjU,n) 


2 

(l+4TT ? ,(*. 2 k 2 ) n+:i/2 


(2.7) 


where k^ is the wavenumber at the location of the jth spectrum 
sample Sj . 

Turning now to derivation of the likelihood equations 
for estimating the spectrum parameters, when Eq. (2.5) is 
substituted into Eq. (3*15) of Ref. 1, one has for the loga- 
of the joint probability density function of the N spectrum 
samples , j=l,2,***,N 


6 


tn[p(S i ,S 2 ,---,S N )] ■ 

« -|N)ln(o 2 L)+)lnF i (Jl,n)+llnF 2 ()l,n) + * • -+*nF N U,n) 

+ [fITTHT + F 2 U,n) + + F N (?,n)J j (2,8) 

where Eq. (3-12) of Ref. 1 also has been used. Equation (2.8) 
is to be maximized with respect to the three parameters 
o 2 L, £, and n. Thus, differentiating Eq. (2.8) with respect 

to o 2 L, x, and n and setting the resulting three expressions 
equal to zero yields after minor manipulations: 


N Sj 

S,f7 

j=i j 


,2 L s I y Lj r 

N AF.(i } n) 


(2.9) 


N 


J-l 


jl UnFjU.n) 


.1 

FjUTnT 


- o 2 L 


= 0 


( 2 . 10 ) 


and 


N 




in * nF j<*-.n) 


S. 

. J. _- 

F,U,n) 

J 


- a z L 


= 0 


( 2 . 11 ) 


Equations (2.9) through (2.11) are three nonlinear simul 
taneous algebraic equations for the spectrum parameters 
a 2 L, l t and n in terms of the periodgram samples 3 2 ,*-*, 

and the functions F^(S.,n) given by Eqs . (2.6) and (2.7) 

for the Bullen transverse and longitudinal spectrum forms 
respectively. However, Eq. (2.9) can be substituted into 
Eqs. (2.10) and (2.11) so that only two nonlinear equations 
need to be solved simultaneously. Let us first define 

G JJt U,n) fcn FjU,n) (2.12) 


7 



9 




(2.13) 


which apply to either the transverse or longitudinal Bullen 
spectral forms as appropriate. Then substituting Eq. (2.g) 
into Eqs. (2.10) and (2.11) and using the definitions 2.12) 
and (2.13)i we obtain instead of Eqs. (2.10) and (2.11), 

E V^fy^br- *E 

j-i L J i=i 




Equations (2.14) and (2.15) are to be solved simulta- 
neously for the two spectrum parameters l and n. These 
solutions l and n are then to be substituted into. Eq. (2.9) 
to yield the parameter o 2 L. Equation (2.4) is then to be 
used to evaluate the scale L, and o 2 then can be evaluated 

by o 2 = o 2 L/L. Thus, Eqs. (2.9), (2.14), and (2.15) are the 
solution to the problem of maximum likelihood estimation of 
the Bullen spectrum parameters o 2 L, £,, and n which yield 
o 2 , L , and n . 

The solution to Eqs. (2.7.4) and (2.15) can be readily 
carried out by numerical means a' ong the lines illustrated 
in Appendix F of Ref. 1. To carry out thes^ solutions, 
the functions G^U,n) and G^ n U,n) are required for the 

Bullen transverse and longitudinal spectra. By differentiat 
ing the natural logarithms of Eqs. (2.6) and (2.7), one 
can readily show for the Bullen transverse spectrum that 


8tt 2 £ 2 k 2 (n+4) 


1— 8tt 2 e, 2 k? (n+1) 


1+ 4 TT 2 1 2 ]^ ( 2 n+ 3 ) + 3 2 IT V 4 k J ( n+ 1 ) 


( 2 . 16 ) 


8 



and 



U,n) 


8ir 2 £ 2 k 2 

*J 

l+8ir 2 * 2 kj(n+l) 


Jlnd+^Tr 2 Jl 2 k 2 ) , 


and for the Bullen longitudinal ep&otrum that 


Gj ^( ^»n) 
and 


1 8tt 2 H 2 k? (n+i) 
— - J ^ 

1 l+lhr 2 Jl 2 kj 


G.„U,n) * - Jtn(l+4TT 2 l z k 2 . ) . 

Jn ,1 


( 2 . 17 ) 


( 2 . 18 ) 


( 2 . 19 ) 
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CONSTRAINED LEAST-SQUARES ESTIMATION OF TURBULENCE 
AUTOCORRELATION FUNCTION PARAMETERS 
First Method 

In this section, the method developed in Sec. 4 of 
Ref. 1 is extended to the case whore the spectrum and auto- 
correlation functions of component w f (t) in Eq. (1.1) are 

better modeled by the Bullen forms rather than the von Karman 
forms used in the cited reference. The parameter estimation 
method developed below is applicable to situations where the 
component w_(t) in Eq. (1.1) is not negligible. 

5 

The autocorrelation function model used here is 

m 

JU) A <t> B U;£,n) + ^Pa^ 1 , 0«US H (3.1) 

1=0 

which replaces Eq. (4.1) of Ref. 1. The quantity cr£ is the 
mean-square value of the ’’fast" component w f ,(t) in Eq. (1.1), 
<pg(S;8,,n) is either the Bullen transverse or longitudinal 

autocorrelation function, as appropriate, with <t> R (0;H,n) * 1, 

and m B 

m 

♦ w U) - a 1 E 1 , (3.2) 

s i=0 


is a polynomial model of the autocorrelation function of the 

"slow” component w (t) in Eq . (1.1) assumed valid within the 

s 

range as in .Sec. 4 of Ref. 1. 

The quantities cj, , £, n, and a., i»0,l,**‘, m are to be 

i X 

evaluated from the autocorrelation function R(C) of a mea- 
sured turbulence record w(t) by minimizing the quantity 



RU) - o*4> B (£;*,n) 


m 



1=0 


2 




(3.3) 
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where the minimization procedure is constrained by the two 
additional relations: 


N 


^°Jn U * n) FjU.n) 


and 





FJT^nT 


y 


(3.4) 


(3.5) 


where periodogram samples and in the wavenumber range 
and k^k^, only are tc be included in the summations 
in Eqs. (3.4) and (3-5) as explained below. 

Ordinarily, the wavenumber spectrum content of the "slow” 
component w (t) in Eq. (1.1) is concentrated in the low 
wavenumber region. Thus, we assume there exists a wavenumber 
k^, such that for wavenumbers kik^,, the *actrum content of 

a turbulence record w(t) is dominated for all k£k f f by con- 
tributions from the "fast” component w^t) only. It is 
advantageous to choose the smallest value of k^, that satisfies 

this condition. The two constraint equations (3*4) and (3*5) 
are to be evaluated using neriodogram samples Sj and taken 

only in the region k£k p f where the wavenumber content of the 

process w f (t) in Eq. (1.1) is everywhere dominant. By 

hypothesis, then, within this region kjk» f , the wavenumber 

spectrum is well represented by the appropriate (transverse 
or longitudinal) Bullen spectral form. Equations (3*4) and 
(3*5) are two of the three likelihood equations available for 
evaluating the three Bullen spectrum parameters. 

To examine why the third likelihood equation (2.14) is 
often not useful in the region k^k^,, we examine the behavior 

of Gj & (£,n) for large values of kj for both the transverse 


n 



<g&***' 



and longitudinal cases. For large values of k*, one can 
readily see from Eq. (2.16) that for the oaee of the Bullen 
transverse spectrum one has 


V**"> - 


SitU’kUn+l) 
J £ 


2 (n+1 ) 
it * 


-8TT 2 Jl 2 k 2 (n+l) 
32 tt ** £** i< ^ (n+1) 


k , 

J 


(3.6) 


and, similarly, from Eq. (2.18) for the case of the Bullen 
longitudinal spectrum one has approximately for large kj 


G JJl (S.,n) 


1 8 tt 2 £ 2 k | (n+^-) 


4irU 2 k 2 


2 (n+1) 


k 

v 


(3-7) 


Thus, for large values of k. for both the transverse and 
;i ongitudinal cases 0^ 0 ( 2. ,n)'- 1 approaches the same constant 

value of -2(n+l)/£. Examination of Figs. 20 and 21 on 

pp. 90 and 91 of Ref. 1 shows that when n * this constant 

value is approximately reached for values of k^l/L. Let 

us now suppose that is approximately equal to 1/L or 

larger. In this case, then, in the region k^k^, where the 

summation in Eq. (2. .VO would be used, we see from Eqs. (3*6) 
and (3.7) that Eq. (2. VO would become, approximately. 


( 2n+l ) 
i 


V' I S 1 l ^ s i 

/^iF^Cn) ‘ N ^ ? T nr7H7 
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or, dividing both sides of Eq. (3.8) by N, 




(3.8) 


which is satisfied for all values of i and n. Hence, when the 
smallest useable value k^ f of kj is in the region where 

1 . 7 ) is satisfied, the likelihood equation 
is satisfied for all values of £ and n and therefore is 
useless for the purpose of determining these values. The 
functions G, (£,n) given by Eqs. (2.17) and (2.19) do not 


or 


jn 

approach constant values with respect to k 


J 


as ^ gets arbi- 


trarily large; hence, the likelihood equation given by Eq. 
(2.15) or (3.4) does not suffer from the above described 
limitation of Eq. (2.1*0. 


To compare the present least-squares methodology with 
that in Sec. *4 of Ref. 1, we see that the Bullen autocorre- 
lation functions $ fl (£;t,n) in Eq. (3-3) have one more parame- 
ter (n) than the von Karman autocorrelation functions 4 >^( 4 ;L) 

in Eq. (*l.3) of Ref. 1. Furthermore, in the present least- 
squares procedure, there is one more equation of constraint, 
Eq. (3- *4), than in the least-squares procedure described in 
Sec. *4 of Ref. 1. Hence, minimization of the quantity ”E" 
given by Eq. (3-3) can be regarded as determining the single 
parameter 7 1, of the "fast” component w^,(t) in Eq. (1.1). 

Similarly, in Sec. 4 of Ref. 1, minimization of "E" given 
by Eq . (4.5) of that reference also determined the single 
parameter o * of the "fast" component Wj,(t). The reader is 

referred to Sec. 4 of Ref. 1 for a more thorough discussion 
of the least-squares procedure. 


Derivation of Algebraic Equations for 

Autocorrelation Function Parameters 

The equation of constraint (3.4) contains two unknown 
parameters i and n. Thus, one can regard that equation as 


13 



determining n as a function o? I — i.e., n»n(i). Hence, 
assuming that this function n=n(M has been determined from 
Eq. (3.4), the second equation of constraint (3*5) can be 
expressed as 


a 


2 


1 

L'[nU> J£N 



(3-9) 


where the subscript f and the overbar have been left off in 
Eq. (3.9), and where the integral scale L has been expressed 
as 

L = L’[n(l)]l , (3.10) 

where from Eq. (2.4), one has 

/tt rf.nU)- 1 -^-] 

L ' [nU;] = — TWTT- (3 - 11) 

where in Eqs. (3.10) and (3.11) n has been expressed as a 
function of l as determined by Eq. (3.4). Hence, Eqs. 

(3.4) and (3.9) together determine a 2 as a function of l — 
i.e., a 2 = a 2 (H). Assuming that o 2 is a monotonic function 
of l, Eqs. (3.4) and (3*9) also determine t as a function of 


l = i(o 2 ) . 


(3.12) 


Thus, assuming that Eq. (3.4) has determined n as a function 
of l, the function <J> B (£;&,n) has been reduced to dependence 

on only a single parameter Si which shall be denoted by 

<t>gU;0 A 4 > B U;t i,nU)] . (3.13) 

Assuming further that i has been determined as a function of 
o 2 = through the additional relationship (3.5), Eq. (3.3) 
then can be expressed as 


14 



where we have used the notation of Eq. (3.14) and the abbre- 
viation o 2 = o 2 . 

The problem, now, is to determine the valuer- of o 2 and 
a^, i=0,l,--*,m that minimize the value of E given by 

Eq. (3*14). This determination is carried out by finding 
value of a 2 and a.,, i=0,i,*-*,m that satisfy 


3E 

3a. 


= 0 , j = 0,1, • • • ,m 


3E 

9a' 


= 0 


(3.15) 


Differentiating Eq. (3.14) with respect to a., j=0,l,*..,m 

7 J 

and with respect to and setting the resulting equations 
equal to zero yields the following set of m+2 nonlinear 
algebraic equations for a , a , •••, a and a 2 : 


0 1 


m 


’H 


m £ 


H • . . 




/ C j 4>B^;Ua 2 )JdC+^a i [C 1+J dC = ^£ J R(£)dS, 

,1 = 0,1, ,m 


i=0 0 


( 3 - 16 ) 


and 




3(j)g:C;^(a 2 )] ( 

— . 


H + 6 R [C;5.(o) 2 [« ;£(a 2 ) ]d£ 


da' 


m 

+ 


9<D n [?,e(a 2 )] 

a — — + Ho 2 )] 

da 2 3 


di 




’H 


r 

-f 


34>b[€U(o 2 )] de 


di 


da 


- + 4> B CC;L(a 2 )] 


R(OdC , (3.17) 
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where the function fc(o 2 ) is determined by Eqs. (3.4) and 
(3.9) as described below. 




% 


Unfortunately, the derivatives of 4>g[€;i(o 2 )] with 


respect to i, turn out to be impractical to evaluate due, in 
part, to the dependence of the parameter n on l that occurs 
in the expression for L'[n(£)] given by Eq. (3.11). Hence, 
the second method described below is recommended instead. 


Second Method 


Let us examine the large k behavior of the Bullen 
transverse and longitudinal spectra, Eqs. (2.1) and (2.2). 
For large k, these spectra are asymptotically equal to 


* BT (k) 


2 T 8tt 2 H 2 k 2 (n+1 ) 
(Hw 2 H 2 k 2 F +3/2 


( 3 . 18 ) 


and 




o 2 L 


(4TrH 2 k 2 ) n+1/2 


(3.19) 


The logarithms of Eqs. (3.18) and (3-19) can be readily 
expressed as 


log$ BT (k)'-log[2(n+l)o 2 L]-(2n+l)log(2TrS,)-(2n+l)log k 


and 


( 3 - 20 ) 


log$> BB (k)''10g(2a 2 L)-(2n+l)log(2iT£)-(2n+l)log (jj.21) 


Hence, when considered as a function of k, Eqs. (3.20) and 
(3.21) are both of the form 


log<J> B (k) 


C - ( 2n+l ) log k. 


( 3 - 22 ) 
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where C is a constant. Thus, when plotted on log- log co- 
ordinates, both the Bullen transverse and longitudinal 
spectral forms are asymptotically linear for large k with 
a slope equal to -(2n+l). In particular, for the von Karman 
transverse and longitudinal spectra, we have n ■ 1/3 which 
yields the well known asymptotic slope of - 5 / 3 * The form of 
Eq. (3.22) suggests the following well-conditioned proce- 
dure for determining the parameter n: 

1. Tabulate and plot the logarithm of the empirically 
determined turbulence velocity spectrum as a function 
of log k. 

2. Identify the region of approximate linearity of the 
data in the large k regime. 

3. Form a least squares fit to the data usirg the 
functional form given by the right-hand side- of 
Eq. (3.22) and solve for C and n. 

4. Retain the value of n thus obtained. 

Once n is evaluated as in steps 1 to 4 above, the Bullen 
spectral forms given by Eqs. (2.1) and (2.2) can be regarded 
as being dependent on only two additional parameters o T and 
L [after substitution of Eq. (2.3) into Eqs. (2.1) and (2.2)]. 
Thus, after evaluation of n as above and substitution of 
Eq. (2.3) into Eqs. (2.1) and (2.2), the procedure for estimat 
ing o* and L given on pp. 92 to 110 of Ref. 1 applies directly 

to the Bullen spectral forms ^^(k) and <I>g^(k) also. The 

only exception is that the von Karman transverse and longitudi 
nal autocorrelation functions and their derivatives given by 
Eqs. (4.48) to (4.52) of Ref. 1 must be replaced by the cor- 
responding Bullen forms which are [5]: for the Bullen trans- 

verse spectra (vertical and lateral components) : 


2 

TTnT 


v\n 


= TTHT ( ^2 ) 


K n (6 ?) - f K n-1 (B ^ 




*b u) ' TOT (¥) 


B5K n (BO-2<n+l)K n _ 1 (Bt) 


for the Bullen longitudinal spectra: 


(3.23) 

(3.24) 
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i 


* . 


♦b«> " rk y (¥) K n (65) 

? B ( ^ ■ * K n-l (8l) > 

where, for both the transverse and longitudinal cases, 
t /rr r(n+^) 

8 • T ‘ - TT n) • • (3 ' 27) 

Equations (3.23) through (3-27) reduce to Eqs . (4.48) through 
(4.52) of Ref. 1 when n is set equal to 1/3 which is the 
value of n for the von Karman case. 


( 3 . 25 ) 

(3.26) 
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